#include<stdio.h>
#include<stdlib.h>
#include<sys/time.h>
#include<math.h>
#include<string.h>
#include<cuda.h>

// Erros
#define ERROR_INVALID_FILE -1
#define ERROR_INVALID_MATRIX -2
#define ERROR_INVALID_VECTOR_SIZE -3

__global__ void divDiag(float *MA, float *MB, float *MADivDiag, float *MBDivDiag, int ordem) {
    int id = blockIdx.x * blockDim.x + threadIdx.x;
	int linha = id / ordem;
    int coluna = id % ordem;
    int idDiagonalVetor = id * ordem + id;
    int idDiagonalMatriz = linha * ordem + linha;

    if(id < ordem) MBDivDiag[id] = MB[id] / MA[idDiagonalVetor];

    if(id < ordem * ordem) {
        if(linha == coluna) MADivDiag[id] = 0;
        else MADivDiag[id] = MA[id] / MA[idDiagonalMatriz];
    }
}

__global__ void calculaXnNext(float *MBDivDiag, float *decrement, float *Xn, int ordem, float *XnNext, float *Xdiff) {
    int id = blockIdx.x * blockDim.x + threadIdx.x;
    if(id < ordem) {
        // Xk+1[i] = B*[i] - sum(A*[i][j] * xk[j])
        XnNext[id] = MBDivDiag[id] - decrement[id];

        // Xdiff[i] = | Xk+1[i] - Xk[i] |
        Xdiff[id] = fabs(XnNext[id] - Xn[id]);
    }
}

// Retorna o maior valor em modulo
float abs_max(float *vector, int size) {
    // Indice do maior elemento
    int abs_max_index = 0;
    // Iterador
    int i;

    // Tamanho invalido
    if(size <= 0) return ERROR_INVALID_VECTOR_SIZE;

    // Percorre o vetor
    for(i = 0; i < size; i++) {
        if(fabs(vector[i]) > fabs(vector[abs_max_index])) abs_max_index = i;
    }

    // Retorno
    return fabs(vector[abs_max_index]);
}

// Rotina principal

int main(void) {
    // Tempo de execucao
    struct timeval tvInitial, tvFinal;
    double tInitial, tFinal;

    // DEBUG
    int debug = 0;

    // Iteradores
    int i, j;
    int numIterations = 0;

    // Variaveis
    int J_ORDER, J_ROW_TEST, J_ITE_MAX;
    float J_ERROR;
    float *MA, *MADivDiag;
    float *MB, *MBDivDiag;
    float *Xn, *XnNext, *Xdiff;
    float *decrement;
    float resultNext;
    char fileName[500], filePath[500];
    FILE *fileInput, *fileResult, *fileDebug;
    size_t vectorSize, matrixSize;
    float errorNext;

    // Variaveis CUDA
    float *dMA, *dMADivDiag;
    float *dMB, *dMBDivDiag;
    float *dXn, *dXnNext, *dXdiff;
    float *dDecrement;
    int blockSize, gridSizeForVector, gridSizeForMatrix;

    // PASSO 1: ABERTURA DO ARQUIVO
    scanf("%s", fileName);

    // Tempo de execucao inicial
    gettimeofday(&tvInitial, NULL);
    tInitial = (double)(tvInitial.tv_sec) + (double)(tvInitial.tv_usec)/ 1000000.00;

	// Abertura do arquivo de entrada
    strcpy(filePath, "\0");
    strcat(filePath, "data/");
    strcat(filePath, fileName);
    fileInput = fopen(filePath, "r");
    if(fileInput == NULL) {
        printf("ERRO! Arquivo invalido!\n");
        return ERROR_INVALID_FILE;
    }

	// Abertura do arquivo de saida
    strcpy(filePath, "\0");
    strcat(filePath, "result_");
    strcat(filePath, fileName);
    fileResult = fopen(filePath, "w");
    if(fileResult == NULL) {
        printf("ERRO! Arquivo de saida invalido (Permissao negada)!\n");
        return ERROR_INVALID_FILE;
    }

	// Abertura do arquivo de debug
	if(debug) {
		strcpy(filePath, "\0");
		strcat(filePath, "debug_");
		strcat(filePath, fileName);
		fileDebug = fopen(filePath, "w");
		if(fileDebug == NULL) {
		    printf("ERRO! Arquivo de saida invalido (Permissao negada)!\n");
		    return ERROR_INVALID_FILE;
		}
	}

    //
    // PASSO 2: LEITURA DOS DADOS
    //

    // Ordem da matriz
    fscanf(fileInput, "%d", &J_ORDER);

    // Fila para ser avaliada
    fscanf(fileInput, "%d", &J_ROW_TEST);

    // Erro permitido
    fscanf(fileInput, "%f", &J_ERROR);

    // Numero maximo de iteracoes
    fscanf(fileInput, "%d", &J_ITE_MAX);

    // TAMANHO DAS MATRIZES
    vectorSize = J_ORDER * sizeof(float);
    matrixSize = J_ORDER * J_ORDER * sizeof(float);

    // ALOCACAO DAS MATRIZES
    Xn = (float*) malloc(vectorSize);
    XnNext = (float*) malloc(vectorSize);
    Xdiff = (float*) malloc(vectorSize);
    MA = (float*) malloc(matrixSize);
    MADivDiag = (float*) malloc(matrixSize);
    MB = (float*) malloc(vectorSize);
    MBDivDiag = (float*) malloc(vectorSize);
    decrement = (float*) malloc(vectorSize);

    // ALOCACAO DAS MATRIZES CUDA
    cudaMalloc((void**) &dXn, vectorSize);
    cudaMalloc((void**) &dXnNext, vectorSize);
    cudaMalloc((void**) &dXdiff, vectorSize);
    cudaMalloc((void**) &dMA, matrixSize);
    cudaMalloc((void**) &dMADivDiag, matrixSize);
    cudaMalloc((void**) &dMB, vectorSize);
    cudaMalloc((void**) &dMBDivDiag, vectorSize);
    cudaMalloc((void**) &dDecrement, vectorSize);

    // LEITURA DAS MATRIZES

    // Matriz A (i x j)
    for (i = 0; i < J_ORDER; i++) {
        for (j = 0; j < J_ORDER; j++) {
            fscanf(fileInput, "%f ", &MA[i * J_ORDER + j]);
        }
    }
    cudaMemcpy(dMA, MA, matrixSize, cudaMemcpyHostToDevice);

    // Matriz B (1 x j)
    for (j = 0; j < J_ORDER; j++) {
        fscanf(fileInput, "%f", &MB[j]);
    }
    cudaMemcpy(dMB, MB, vectorSize, cudaMemcpyHostToDevice);

    // DEBUG: Imprime os dados
    if(debug) {
        fprintf(fileDebug, "DADOS DO PROBLEMA\n");
        fprintf(fileDebug, "Ordem da matriz: %d\n", J_ORDER);
        fprintf(fileDebug, "Fila para ser avaliada: %d\n", J_ROW_TEST);
        fprintf(fileDebug, "Erro permitido: %f\n", J_ERROR);
        fprintf(fileDebug, "Numero maximo de iteracoes: %d\n\n", J_ITE_MAX);
        fprintf(fileDebug, "MATRIZ A\n");
        for (i = 0; i < J_ORDER; i++) {
            for (j = 0; j < J_ORDER; j++) {
                fprintf(fileDebug, "%f ", MA[i * J_ORDER + j]);
            }
            fprintf(fileDebug, "\n");
        }
        fprintf(fileDebug, "\n");
        fprintf(fileDebug, "MATRIZ B\n");
        for (j = 0; j < J_ORDER; j++) {
            fprintf(fileDebug, "%f ", MB[j]);
        }
        fprintf(fileDebug, "\n\n");
    }

    // DEFINICAO DOS BLOCOS DE THREADS
    blockSize = 500;
    gridSizeForVector = (int) ceil((float) J_ORDER / blockSize);
    gridSizeForMatrix = (int) ceil((float) (J_ORDER * J_ORDER)/ blockSize);

    //
    // PASSO 3: EXECUCAO DO ALGORITMO
    //

    // Encontra as matrizes MADivDiag e MBDivDiag
    divDiag<<<gridSizeForMatrix, blockSize>>>(dMA, dMB, dMADivDiag, dMBDivDiag, J_ORDER);
    cudaMemcpy(MADivDiag, dMADivDiag, matrixSize, cudaMemcpyDeviceToHost);
    cudaMemcpy(MBDivDiag, dMBDivDiag, vectorSize, cudaMemcpyDeviceToHost);

    // DEBUG: Imprime os dados
    if(debug) {
        fprintf(fileDebug, "MATRIZ A*\n");
        for (i = 0; i < J_ORDER; i++) {
            for (j = 0; j < J_ORDER; j++) {
                fprintf(fileDebug, "%f ", MADivDiag[i * J_ORDER + j]);
            }
            fprintf(fileDebug, "\n");
        }
        fprintf(fileDebug, "\n");
        fprintf(fileDebug, "MATRIZ B*\n");
        for (j = 0; j < J_ORDER; j++) {
            fprintf(fileDebug, "%f ", MBDivDiag[j]);
        }
        fprintf(fileDebug, "\n\n");
    }

    // Solucao inicial
    for(i = 0; i < J_ORDER; i++) {
        Xn[i] = MBDivDiag[i];
    }


    // Metodo iterativo
    do {
        // Transpoe Xn para o CUDA
        cudaMemcpy(dXn, Xn, vectorSize, cudaMemcpyHostToDevice);

        // Calcula a cumulativa
        for(i = 0; i < J_ORDER; i++) {
            decrement[i] = 0.0f;
            for(j = 0; j < J_ORDER; j++) {
                decrement[i] += MADivDiag[i * J_ORDER + j] * Xn[j];
            }
        }
        cudaMemcpy(dDecrement, decrement, vectorSize, cudaMemcpyHostToDevice);

        // Calcula XnNext e Xdiff
        calculaXnNext<<<gridSizeForVector, blockSize>>>(dMBDivDiag, dDecrement, dXn, J_ORDER, dXnNext, dXdiff);
        cudaMemcpy(XnNext, dXnNext, vectorSize, cudaMemcpyDeviceToHost);
        cudaMemcpy(Xdiff, dXdiff, vectorSize, cudaMemcpyDeviceToHost);

        // Calcula o valor aproximado resultado do processo
        resultNext = 0.0f;
        for(j = 0; j < J_ORDER; j++) {
            resultNext += MA[J_ROW_TEST * J_ORDER + j] * XnNext[j];
        }

        // Calcula o erro
        errorNext = abs_max(Xdiff, J_ORDER) / abs_max(XnNext, J_ORDER);
		
		// Atualiza a solução
		for(i = 0; i < J_ORDER; i++) {
			Xn[i] = XnNext[i];
		}

        // Atualiza o numero de iteracoes
        numIterations++;
    }
    while(errorNext >= J_ERROR && numIterations < J_ITE_MAX);

    // Imprime os resultados do processo no arquivo de resultado
    fprintf(fileResult, "Iterations: %d\n", numIterations - 1);
    fprintf(fileResult, "RowTest: %d => [%f] =? %f\n", J_ROW_TEST, resultNext, MB[J_ROW_TEST]);

    // Tempo de execucao final
    gettimeofday(&tvFinal, NULL);
    tFinal = (double)(tvFinal.tv_sec) + (double)(tvFinal.tv_usec)/ 1000000.00;

    // Imprime o tempo de execucao no arquivo de resultado
    fprintf(fileResult, "Tempo de execucao: %.2lfs\n", (tFinal - tInitial));

	// Imprime os resultados do processo na tela
    printf("Iterations: %d\n", numIterations - 1);
    printf("RowTest: %d => [%f] =? %f\n", J_ROW_TEST, resultNext, MB[J_ROW_TEST]);

	// Encerra os arquivos de entrada e saida
	fclose(fileInput);
	fclose(fileResult);
	if(debug) fclose(fileDebug);

	// DESALOCA OS PONTEIROS
	free(Xn);
	free(XnNext);
	free(Xdiff);
	free(MA);
	free(MADivDiag);
	free(MB);
	free(MBDivDiag);

    // DESALOCA OS PONTEIROS CUDA
    cudaFree(dXn);
    cudaFree(dXnNext);
    cudaFree(dXdiff);
    cudaFree(dMA);
    cudaFree(dMADivDiag);
    cudaFree(dMB);
    cudaFree(dMBDivDiag);

   // Retorno
    return 0;
}
